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We present a general method based on nonlinear response theory to obtain effective interactions 
between ions in an interacting electron gas, which can also be applied to other systems where an 
adiabatic separation of time scales is possible. Nonlinear contributions to the interatomic potential 
are expressed in terms of physically meaningful quantities, giving insight into the physical properties 
of the system. The method is applied to various test cases and is found to improve the standard 
linear and quadratic response approaches. It also reduces the discrepancies previously observed 
between perturbation theory and density functional theory results for proton-proton pair potentials 
in metallic environments. 



I. INTRODUCTION AND BACKGROUND 

Perturbation and response theories are ubiquitous in 
physics. Here we will discuss response theories as ap- 
plied to many-component systems, where the degrees of 
freedom associated with one component can be traced 
out to yield, via an adiabatic separation of time scales, 
effective interactions between the remaining compo- 
nents. These effective interactions have had their share 
of success in describing simple metals 1,2 , metals with 
magnetic impurities [i.e. the Rudcrman-Kittcl-Kasuya- 
Yosida (RKKY) interaction 3 ], or, in the classical do- 
main, colloidal suspensions [e.g. the Derjaguin-Landau- 
Verwey-Overbeek (DLVO) interaction between charged 
colloidal particles 4,5 or the depletion interaction 6,7 ]. Lin- 
ear response in particular has been used extensively to 
obtain pair potentials; its simplicity makes it a very valu- 
able qualitative tool, and in many contexts it also pro- 
vides sufficient quantitative accuracy. 

Higher-order contributions have been used less often, 
partly because calculational complexity increases rapidly 
with the order of perturbation considered. The first 
derivation of the static quadratic response function for 
the ground-state noninteracting electron gas originates 
with Lloyd and Sholl 8 . It was later rederived by Brov- 
man, Kagan et al, who wrote a series of articles ap- 
plying nonlinear perturbation to metallic systems (see, 
for example, Brovman, Kagan, and Kholas 9 , and ref- 
erences therein). A third derivation of the static re- 
sponse function was provided even later by Milchev and 
Pickenhcim 10 . A derivation of the dynamical quadratic 
response function can also be found in Pitarke et al 11 . 
With the advent of density functional theory (DFT) 
and modern computers, response theory has lost some 
of its quantitative appeal, but as noted it is still being 
used for a qualitative understanding of general trends in 
metals 11-13 . 

Dynamic compression experiments have achieved met- 
allization of hydrogen at high temperature 14 , and there is 
an expectation that low-temperature metallization might 
also be achieved in the relatively near future 15 . Hydro- 
gen has always been difficult to treat by perturbation 
approaches because the interaction of a proton with the 
electron gas is not weakened by a repulsive core, as is 



the case for many elements 16 . It is not even obvious that 
a perturbation approach converges. At high pressures, 
though, the enhanced kinetic energy of the electron gas 
in hydrogen leads to a better convergence of response 
theory. Response-based pair potentials were in fact used 
to describe with some success the pressure dependence 
of the vibron in the hydrogen solid 17 . On the other 
hand, the application of standard DFT methods to hy- 
drogen at very high densities is challenging because of the 
strongly inhomogeneous, cusplike electronic density sur- 
rounding the proton and the possible failure of pseudopo- 
tcntials designed for a different density range. The pair 
potentials obtained from quadratic perturbation theory 
for dense hydrogen exhibit discrepancies when compared 
with ab initio DFT results 18 for Wigner-Seitz radii as low 
as 1.3. The Wigner-Seitz radius r s is related to the un- 
perturbed electronic density po by r s a = [3/(47r/9 )] 1 ^ 3 , 
where a — h 2 /m e e 2 is the Bohr radius. 

In order to use the perturbation method as a reli- 
able, analytic alternative approach to DFT in this den- 
sity range, higher-order terms are needed. Going beyond 
quadratic response using conventional methods is a sub- 
stantially more difficult task. Accordingly, we derive a 
simple identity that allows us to obtain effective pair po- 
tentials beyond quadratic response, which amounts to 
carrying out a partial sum of perturbation terms. The 
result is not only very intuitive, expressing the pair po- 
tential in terms of physically meaningful expressions, but 
also very general, since it applies in any dimension and for 
both classical and quantum systems. Moreover, it does 
not require explicit knowledge of the nonlinear response 
functions. It should be noted that if nonlinear effects 
are to be fully taken into account, many-body interac- 
tions should also be included, either directly or through 
effective-medium approaches 33 . This article focuses on 
obtaining pair potentials that can be used as a starting 
point for cither approach. 

In Section II, we outline the response theory formalism 
and derive the simplest version of our result, equation 12. 
In Section III (and Appendix I) a generalization is de- 
rived that applies to homogeneous, interacting systems. 
The intuitive nature of this result permits us to explain 
part of the discrepancies observed between the results of 
Nagao et al 17 and Bonev and Ashcroft 18 , and to improve 
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upon quadratic response. In Section IV we apply our re- 
sults to obtain pair potentials in various systems, includ- 
ing metallic hydrogen. Generalizations to many-center 
interactions and magnetic perturbations are discussed in 
Appendix II. 



II. EFFECTIVE INTERACTIONS AND 
RESPONSE THEORY 

We start with a canonical system composed of a mix- 
ture of at least two different types of particles, in con- 
tact with a heat bath establishing a temperature T. The 
ground-state properties of the system can be studied from 
the limit T —> 0. 

The particles will be divided in two subsets, labeled 
L (for Light, typically electrons) and M (for Massive, 
typically ions). It is assumed that all L particles are 
identical. On the other hand, the set M can contain 
various types of particles. The time scales associated 
with the L particles will be assumed to be much shorter 
than those associated with the M particles, allowing the 
use of an adiabatic separation of time scales. In real 
systems, this separation is not exact. The assumption 
of exact separation of time scales, which will be made 
throughout this article, is referred to as the adiabatic or 
Born-Oppcnheimer approximation. We will not discuss 
the accuracy of this approximation here and refer the 
reader to the existing literature (see, e.g., Ziman 19 ). 

The adiabatic approximation allows the treatment of 
the degrees of freedom associated with the light particles 
L in an average way, leading to a much simpler effective 
Hamiltonian for the remaining M particles. 

The Hamiltonian of the initial system has the form 



i 3 

+ V iM ({r},{R}). 



The notation {r} designates the set of all coordinates 
{^i}i=i,...,N L i and Pi,i"i, and m are the momenta, posi- 
tions and mass of L particles, while P i; Rj, and are 
the corresponding quantities for M particles. We imag- 
ine the system to be confined to a macroscopic volume 
ft. 

In many relevant physical systems the interactions be- 
tween the particles are largely pairwise. An important 
example of such a system is a metal, where the L parti- 
cles could be the valence electrons and the M particles 
the ionic cores. For pairwise systems we can write 



V L ({r}) = J2^\r i )+^v L (r i ,r j ), 

i i=jLj 

V M ({R}) = ]T 0#(RO + J2 »« (Ri.B-i). (2 ) 

i i=jLj 

V LM ({r},{K}) = J2vf M (r i ,K j ), 

ij 

where 4>l^ an d </>2r denote external, one-particle poten- 
tials, v L is the pair potential between L particles, vfj 
is the pair potential between M particles i and j, and 
v^ M is the pair potential between the M particle j and 
the L particles. In the case of long-range, Coulombic 
interactions, care should be taken in the definition of 
the pair potentials to ensure that both the unperturbed 
and the perturbing systems are thermodynamically well 
defined 20 . 

We will restrict our attention to such pairwise systems, 
even though our main results require only the slightly 
weaker assumptions that we can write the interaction 
between L and M particles as 

V LM ({r} 1 {R}) = J2v l LM (r t ,{R}). 

i 

This is to say that the joint effect of all the M particles, 
if they were held fixed, would amount to an additional 
one-body potential for the L particles. 

A. Effective interactions and the adiabatic 
separation of time scales 

As mentioned, we want to find an effective Hamilto- 
nian for the M particles by tracing out the degrees of 
freedom associated with the L particles. Since volume 
(ft) and temperature (T) are specified, we begin with 
the relevant Hclmholtz free energy F, 

F = -fcTlnTr {LM} e- /3H '. (3) 

We then trace over the degrees of freedom associated 
with L at fixed configuration of M, assuming an adiabatic 
separation of time-scales, 

F = -fcTln (Tr M e^ TM+yM ^r L{M) e-^ TL+vL+vLM ^ , 

where T^um) means the trace over states of particles 
of type L for a fixed configuration of particles M, and 
rpL,M j g j. ne k me ti c energy associated with the L and M 
particles, respectively. The free energy of system L for a 
fixed configuration of M is simply 

F L ({R},T) = -fcTln (Tr i(M) e-^ L +^" f >) . (4) 
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The total free energy can then be written as 

F = -kTlnTv M e- f3H M f 7 (5) 
where the effective Hamiltonian has the form 

Heff = E + ^ M ({ R i) + v eff ({n}) 

i 3 

and V eff ({R}) = F L ({R},T) is the desired effective in- 
teraction between particles of type M. It is clearly state 
dependent, since it depends on temperature. It also de- 
pends on the volume and on the properties of the L 
particles, including their mass, their number, and the 
form of their interactions. Note that even if the ini- 
tial Hamiltonian contains only pairwise interactions, the 
effective Hamiltonian will typically contain many-center 
interactions as well as volume-dependent but structure- 
independent terms. Here we will focus on the determi- 
nation of pair interactions. 



The functions are i by definition, the response 

functions of the unperturbed system (which at this point 
is not necessarily uniform) and carry all the information 
about this system, including temperature. We have as- 
sumed a large but finite volume ft and a dense but dis- 
crete distribution of wavevectors k, since our main in- 
terest is the electron gas. Continuous systems can be 
recovered using the usual prescription 20 [in three dimen- 
sions, it reads ~~ * rfk/ (2tt) 3 ] . Even though the 
response functions depend on temperature, the methods 
we present here do not involve this temperature depen- 
dence explicitly In order to simplify the notation, we 
will not explicitly keep track of the temperature in the 
following. Unless otherwise specified, xW(k, k') will sim- 
ply stand for x (1) ( k , k', T), etc. 

Equation (6) can be used together with the coupling- 
constant integration method 21 ' 22 to obtain the variation 



In order to obtain the effective interaction we need to 
calculate the free energy (4) of system L while the M 
particles are held fixed. To do this we treat the poten- 
tial Vf M (r, {R}) as an external one-body perturbation 
to the system composed of particles of type L. In the fol- 
lowing Vj LM (r, {R}) will therefore be simply written as 
Vext( r )> keeping in mind the dependence of V ext on the 
Ri. We then proceed to a functional expansion of the 
free energy in orders of V ext (k) w 

V ext (k) = f dre- ik - r V ext (r). 
Jn 



B. Response theory, and the coupling constant 
integration method 

The change in density induced in a system by a per- 
turbing external potential V ext (k) takes the form 



(6) 
I 

in the Helmholtz free energy arising from the perturba- 
tion V ext : namely, 



AF= f 1 dX{V ext ) x , (7) 
Jo 

where (-) A is the statistical average with respect to the 
states of the Hamiltonian H^ xt , which is obtained by re- 
placing V ex t by XV ex t in H ext . In the thermodynamic 
integration scheme 23 , the integrand in equation (7) is de- 
termined by numerical simulation. In the perturbation 
approach, it is instead expanded in powers of the external 
perturbation, yielding 



Sp(k,T)=p(k,T)-p^\k,T) 

= ^x«(-k,k',T)v ext (k') + ^ ]T X (2) (-k,k',k",T)y ext (k')^t(k") + • • • 

k' k',k" 



k,n 



P ^(k)V ext {-k) 
n+l 



= E 



- (n + l)fi« 



J2 X (n) (ki , k„ +1 ) V ext (kx ) . ..V^t (k„+i ) , 



(8) 



n=0 



ki,...,k„ 



where p^ is the part of (6) that is of order n in V ext . The of the unperturbed system, x^°^(k) = p^°\— k)/f2. 
"zeroth order" response function is related to the density 
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Note that terms of n order in the external potential 
in expression (6) yield terms of order n + 1 in (8). In the 
following, " n th order response" refers to terms of order n 
in (6), unless otherwise specified. 

In order to study pair potentials, let us now as- 
sume that the perturbation originates with two external 
sources, located at positions R and R&, so that 

Vext(k) = K(k)e ik ' R * + V b (k)e^ R \ 

The free energy now depends on R a and R& (and would 
depend only on R a — R h if we had restricted ourselves 
to initially homogeneous systems). This yields an in- 
duced effective pair potential 0*(R a , R&) between the two 



where Sp\ m (— k, R^) is the density induced, at linear or- 
der, by particle i located at position R^. 

The first two terms in this expression correspond to 
the Coulombic energy of the system at the level of linear 
response. The third term, which arises from the coupling- 
constant integration, incorporates the variation in kinetic 
energy and entropy caused by the perturbation. 

Instead of keeping terms linear in V ext , we can choose 
to keep all terms that are linear in V a , for all orders in 
Vb- In this case we find, using equation (6), that the pair 
potential can be written as 

</> 4 (R a ,R fc ) ~ 0g Pb (R a , R b ) 

= ^E^( k )^(- k ' R ^ kRa - 

Here, dpb(k, R&) is the total density that would be in- 
duced if the perturbation potential was caused by source 
b alone, i.e., V ext (k) = Vb(k)e lk ' Rb . Notice that the terms 
arising from the interaction between 5p a and Vb and the 
change in kinetic energy and entropy arising from the per- 
turbation (taken into account by the coupling-constant 
integration) cancel each other out exactly in this case. 
One can also obtain a different estimate to the potential 
by inverting the roles of a and b in the previous discus- 
sion, leading to <p l SPa . 

This asymmetric approach will be referred to as the 
successive perturbation method, or SPM. Indeed, one in- 
terpretation of this result is that instead of perturbing 
the initial system with a and b simultaneously, we ini- 



sources, defined as the sum of all terms in (8) that de- 
pend on both V a and V b - Here the standard approach is 
to keep only terms up to a given order in V ext - Linear 
response yields the induced pair potential 

4„(R ,R fe ) = i£x (1) (k,k') 

k,k' (9) 

x K(k)y 6 (k')e i(k ' Ra+k '' Rb) - 

In order to emphasize the role of the coupling-constant 
integration and to hint at an upcoming result [equation 
(12)], this can be written as 



(10) 



I 

tially add only 6, determine the properties of the inter- 
mediate system, and then, subsequently, add particle a. 
Since our pair potentials result from the calculation of a 
free energy, the Gibbs-Bogoliubov inequality applies and 
equation (11) is related to a rigorous bound on the pair 
potential: namely, 

4 (Ra, R&) < 0sp 6 (R a , R 6 ) + (K)o - AF a , 

where (V a )o is the statistical average of the operator 
V a (k)e 4k Ra with respect to states of the initial, unper- 
turbed system, and AF a is the change in free energy 
induced by adding a to the unperturbed system, which 
can be obtained by setting Vb=0 in (8). 

Equation (11) is the first term in the expansion of the 
energy if b is treated exactly and a is treated as a pertur- 
bation. Note that higher-order terms in the expansion 
would simply involve higher-order functional derivatives 
of the free energy of the intermediate system including b, 
which are related to density-density correlation functions 
of that system. 

The SPM is especially useful if V a is weak and Vb is 
strong 24 , but is not as useful when both sources of per- 
turbation are strong and require nonlinear treatment. On 
the other hand, we can easily obtain from this relation 
another expression for the induced pair potential that 
incorporates all contributions that arc linear in either 
source. This expression therefore includes all contribu- 
tions to the pair potential up to cubic order in V a and 
Vb. It reads 



^ m (R a ,R h ) =- ]T (K(k)C(-k,R i )e ,k ' R ' + V b (k)6p l : n (-k,R a y^) 

k 

-JE X (1) (ki,k 2 )K(k 1 )H(k 2 ) e ^ kl - R » +k2 - R ^, 
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4 S p(Ra,R 6 ) =^(K(k)<5p b (-k,R fc ) e 4k ' R » + V b (k)Sp a (-k,R a )e^) 

I (12) 
- X (1) (ki,k 2 )K(k 1 )H(k 2 )e l ( kl - R "+ k - R ^. 



ki,k 2 



This result, which emerges from what we will refer 
to as the symmetrized successive perturbation method 
(SSPM), can also be seen to be a natural generalization 
of linear response (10). Equation (12) is our main result 
for noninteracting systems. It also applies to interacting 
systems, but in that case it can be improved upon. We 
will do this in Section III. 

Equation (12), despite its simplicity, has many inter- 
esting features: 

(a) It is very general; it can be used in any dimension, 
for classical and quantum, homogeneous and inhomoge- 
neous, noninteracting and interacting systems (although, 
as mentioned, it can be improved for interacting systems; 
see Section III), (b) It is intuitive: the first two terms 
are the interaction of the potential energy associated with 
each perturbation with the density induced by the other. 
The last term, which is equal to the negative of the lin- 
ear response potential, accounts for the change in kinetic 
energy and entropy of the perturbed system and the con- 
tributions to the density that are not additive in V a and 
Vb- (c) It includes all contributing terms up to quadratic 
response [yielding cubic terms in equation (8)], plus 8 out 
of 14 contributing terms of cubic response: it includes 
more terms than quadratic response, (d) It does not re- 
quire the explicit knowledge of the second-order response 
function, but only that of the linear response function of 
the initial, unperturbed system, (e) It expresses the pair 
potentials in terms of quantities that are simple, sym- 
metric, and can in principle be measured. Finally (f), 
since it takes the effect of a single, isolated perturbation 
as an external input, equation (12) allows us to treat the 
effect of stronger, localized perturbations which are often 
difficult to treat with purely perturbative methods. 

We use this result in two test cases in section IV. The 
first example is the effective interaction between particles 
perturbing, via delta-function potentials, a noninteract- 
ing quantum one-dimensional electron gas. The second 
is a version of the classical Asakura-Oosawa model 6,7 of 
the depletion interaction, with finite square wells replac- 
ing hard-sphere potentials. We then apply it to the more 
realistic calculation of the pair potential between protons 
in a metallic environment. 

Before we do this, we use the intuitive form of equation 
(12) to suggest the existence of a higher-order correction 
that applies to interacting systems, which will be derived 
in Appendix I. 

III. CORRECTIONS SPECIFIC TO 
INTERACTING SYSTEMS 

As mentioned, equation (12) is valid for the interact- 
ing electron gas as well as for the noninteracting one: 



the interactions simply modify the forms of the ■ But 
upon further inspection of this equation, one might won- 
der about the absence of coupling between the induced 
densities themselves. The part of this interaction that is 
linear in V a or Vb is included in (12), but the part that is 
nonlinear in both V a and Vb is not. It turns out that some 
of these contributions can also be expressed intuitively in 
terms of Sp a and Spb- 

For homogeneous interacting systems we obtain, by 
summing up higher-order terms that correspond to re- 
ducible diagrams (see Appendix I) , an extra term of the 
form 



«Ra&) = ^^r(k)^W^f L (-k) e 4k - R -, (13) 
k 

where H ab is the distance between a and b, Sp^ L — 
S p a — fipa \ and v is the effective interaction between 
the induced densities. It reads 

v(k)=e(k)c^\k), (14) 

where 

e(k) = l-c$\k) X0 (k) 

involves the noninteracting response function Xo^- ^ n 
classical statistical mechanics Cq (fc) is the Ornstein- 
Zcrnikc function for the unperturbed system, multiplied 
by k B T. In the case of quantum mechanical particles of 
charge e, we have 20 

c£\k)=v c (k) + n 1 (k), 

where v c (k) — Aire 2 /k 2 is the Coulomb interaction and 
pi is the first functional derivative with respect to density 
of the exchange-correlation potential. It is related to the 
local-field correction G by pi(k) = —G(k)v c (k). 

With some reorganization, we can now write the to- 
tal induced pair potential in a remarkably simple form: 
namely, 



G 



= < k ) V a (k)8 Pb (-k) + 5p a (k)V b (-k) + 5p a (k)c£\k)Sp b (k) - V a (k) x ( o>(k)V b (-k) 



(15) 



In real space, 



'(R-ab) 



rE 1 



! (k)e 



tk-R. 



This result brings up an important point. The energy 
associated with the interaction between potential V b and 
the n th order contribution to the density induced by a, 

which we write as Sp a n \ arises at order n+l in the energy 
expansion. On the other hand, the "electron-electron" in- 
teraction between Sp a n ^ and Sp^ appears only at order 
2n in the energy. Except for the special case of linear re- 
sponse, the electron-electron interaction terms originate 
with a higher order in response than the corresponding 
electron-perturbation terms. Therefore termination of 
the series (8) at any order beyond linear response will 
typically result in pair potentials that are overly attrac- 
tive. 

Indeed, if we compare the effective pair potentials 
for hydrogen atoms in jcllium obtained from quadratic 
response 17 to those obtained from ab initio methods 18 , 
we observe exactly such a discrepancy. We calculate in 
section IV C the lowest-order contribution arising from 
(13), and find that it indeed improves the agreement be- 
tween ab initio methods and response theory. 

We finally draw the reader's attention to the similar- 
ity (for protons in an electron gas) between the ab initio 
equation (15) and the variational Hcitlcr-London evalu- 
ation of the isolated hydrogen molecule energy. This will 
be further discussed in section IV C. 



IV. EXAMPLES OF APPLICATIONS 



A. Delta functions in a noninteracting electron gas 



Sp Un {r) 



2m 



J2x W (k,-ky kr (17) 



and is the linear response function of the one- 

dimensional noninteracting electron gas 3,25 : 



fei + 2k F 



k\ — 2k F 



5fei,-fc 2 - 



The Fermi wave-number kp is related to the unper- 
turbed linear density by kF = irpo- 

Converting the sum into an integral, we find (see 
Kittel 3 , and also Yafet 25 and Giuliani et al 26 ) 



Sp Hn (r) = -Si(2k F r) - " 

7T I 



(18) 



where Si is the sine integral function 27 . The nonlinear in- 
duced density Sp a (r) can also be calculated exactly as the 
sum of the bound and scattering state density: namely, 



Spa = Spbound + Sp s 



with 



(19) 

Here 6 is the Heaviside step function. 
The scattering state density is, in the limit of large 
L 2e 



As a simple instructive example, we first consider a 
one-dimensional noninteracting electron gas of unper- 
turbed linear density po confined to a large length L, 
with periodic boundary conditions. The external pertur- 
bations have the form 



V «(r) = ^5(r) (a = a,b). 

From (12) we obtain immediately 



bssp(Rab) - ^(2S Pa (R ab )-Sp Un (R ah )). 



Here 



(16) 



Jo 



kF (2uksin{2k\r\) u 2 cos(2fc|r|) 



Ak 2 



Ak 2 + u 2 



(20) 

Factors of 2 are included for spin degeneracy. This inte- 
gral can be evaluated in terms of the exponential integral 
Ei using, for example, relations 5.1.41 and 5.1.42 from 
Abramowitz and Stegun 27 . We find 



e u\r\ u 



(3m [Ei ((u + 2ik f ) \r\)] + tt6(-u)) . 



(21) 

The total density induced by a single delta function 
potential in a noninteracting electron gas therefore takes 
the very simple form 
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ue 



,\r\ 



-3m [E x ((it + 2%) |r| 



(22) 



Note that for kf —* we find that <5p = Spbound, as 
expected. The appearance of a bound state at u = 
corresponds to the branch cut of Ei along the negative 
real axis. 

If we use (18) and (22) in (16), we find an expression 
for the pair potential as a function of R a b which reads 



4>SSp(Rab) 



2, ,2 



h u 



27TTO 



„uR a 



(3m [E 1 ((u + 2%) J R ob )D 



(si(2k F R ab ) - |) 



(23) 



Alternatively, the exact pair potential can be calcu- 
lated by solving the Schrodingcr's equation directly for 
different values of R a b- The result is expressed as a 
sum over the eigenvalues of the Hamiltonian, which are 
obtained by numerically solving a set of transcendental 
equations. 

We compare in Figure 1 these numerical pair poten- 
tials with the analytical ones obtained from the SPM and 
SSPM [equation (23)], and those obtained from linear 
response, for various interaction strengths, including at- 
tractive and repulsive cases. We observe that the SSPM 
improves upon linear response and the SPM, especially 
quantitatively, at low u, but also qualitatively (at higher 
u.) 

As a side remark, notice that the system has no bound 
state for u > 0, one even bound state for u < 0, and 
one extra odd bound state when R ab u < —2. It has 
been suggested 28 that the appearance of a bound state 
might cause the failure of response theory, since such 
states are qualitatively different from the initial, unper- 
turbed free-electron system. We see that this is not the 
case in this particular example and that the pair poten- 
tials can be accurately described by response theory de- 
spite the presence of a bound state. This is consistent 
with the observation that no discontinuity in the density 
of charge occurs upon the formation of a localized, bound 
state (see, e.g., Galindo et at 29 ). 



B. Classical depletion interaction 

The depletion interaction (or entropic attraction) plays 
a major role in classical colloidal systems, where it is 
used to tune the interaction between colloidal particles 
immersed in a solvent. The addition of polymers to the 
solvent indeed causes an effective attraction between the 
colloidal particles which can be adjusted by modifying 
the polymer concentration. 

The Asakura-Oosawa model 6,7 describes this effect by 
representing colloidal particles as hard spheres of radius 
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FIG. 1: (color online) Pair potentials for delta function po- 
tentials in one dimensional noninteracting electron gas: exact, 
linear response, and the SPM and SSPM results. The Fermi 
wavevector is set to kf = 0.785/ao. The residuals (difference 
between estimated and exact pair potentials) are shown in 
the inset. Note the difference in behavior of the SSPM near 
the origin between attractive (first three figures) and repulsive 
(last figure) delta function potentials. 
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D and (folded) polymers as hard spheres of radius 5. The 
interactions between the polymers are neglected. The 
polymers are therefore treated as an ideal gas which is ex- 
cluded from spheres of radius D + S surrounding each col- 
loid. Here we will replace the hard-sphere potential with 
a finite repulsion of magnitude Vq. Within this model the 
effective attractive potential between the colloids can be 
calculated exactly and compared to the response theory 
results, providing a useful benchmark for perturbative 
approaches in the classical regime. 

Suppose we have only two colloidal particles (at posi- 
tions R a ,t>) in a bath of polymers of unperturbed density 
po = N/£l. Since the polymers are taken not to interact, 
their density is given by 



P(r) 



Ne -0V e (r) 



where V e (r) = V a (r) + V b (r) is the total potential and 
_ / V Q if |r - R a , 6 | < D + S, 



(24) 



V a , 



otherwise. 



Since V e (r) is zero except within a bounded region, 
p(r) simplifies, in the thermodynamic limit, to 



(25) 



Now we can use the coupling-constant integration 
method to obtain the Hclmholtz free energy: 



F-F Q =J^ d\J dvV ext (v) Px (v) 



= -A (R a ,R b )^(e- 2 ^-l) (26) 

- Ar (R a ,n b )^(e-e v °-i), 

where A Q is the overlap volume between the two spheres 
of radius D + S centered at R ,b, and A r = 8tt(D + 
5) 3 /3 — 2A is the non-overlapping volume of the spheres. 
All the R ai f, dependence is contained in A (R a ,Rb) = 
A (R a — Rfc). Accordingly we write the pair potential as 



l (R. 



a R b , v ) = _ p°M*«-**) {e -m 1)2 



(27) 

Notice the obvious limit V — > oo, yielding the familiar 
Asakura-Oosawa result 6,7 



^(R a -R b ,oo) = - ^ o(R ^- Rb) . (28) 

This problem can also be treated at various orders of 
response theory, using, e.g., x^C 4 , k') = — /3po^k,-k'- 



The general form for the pair potential, which can be 
obtained by carrying the perturbation to any order, can 
be expressed as 



bl(R a -R b ,Vo) = -f(pV ) 



p A (R a - Rj) 



Only the position-independent function / is modified 
in the various approximations. The exact (nonperturba- 
tive) form for the function / is 

f{x) = (e-*-lf. 

We evaluate the performance of the various orders in 
response and of the SPM and SSPM by comparing the 
estimates obtained from each method to the exact result 
for f(x). Since we know the exact result analytically, we 
can obtain the n th order response estimate to /, and 
hence to the pair potential, by a simple Taylor expansion 
of f(x), keeping terms up to (n + l)*' 1 order in x = (3Vo, 
i.e., 



/(*) = x 2 - x 3 + -x 4 + ■ ■ ■ 

Expansions up to cubic order are shown on Figure 2, 
together with the SPM result [equation (11)] 



fsPM(x) = (1 - e~ x ) ^^ 2 -y + y + -- - 

and the SSPM result [equation (12)] 



Isspm{x) = 2 ( x{l - e x ) - — 



x A - x> + — + 



Negative values of x are included in Figure 2 for il- 
lustrative purposes. One can again observe that because 
of the additional terms it includes, the SSPM is more 
accurate than linear and quadratic response for all V a 
and more accurate than the SPM for x < 1.59. Cubic 
response, on the other hand, is closer to the exact result 
than the SSPM for -1.93 < x < 0.78, which is expected 
since the SSPM does not include all cubic terms. For 
larger \x\, though, the higher-order terms play a more 
important role and the SSPM is more accurate than cu- 
bic response. 

This example is also instructive in that it allows us to 
study directly the convergence of response theory. Note 
that for Vo < the convergence is monotonous, while for 
Vo > there is a more complex, alternating approach to 
the exact result. This can be traced back to the fact that 
changing the overall sign of the perturbing potential re- 
sults in changing the sign of odd orders in response, with- 
out affecting the even orders. Thus, if a potential exhibits 
monotonous convergence, its additive inverse exhibits al- 
ternating convergence. We might therefore expect, for 
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FIG. 2: (color online) (upper figure) Dependence of the effec- 
tive pair potential between two repulsive colloidal particles on 
the interaction strength parameter x = /3Vo- exact result ver- 
sus the SPM and SSPM, and perturbations at various orders, 
(lower figure) Relative errors \ f n (x) — f(x)\/f(x). 

example, that effective interactions with protons and an- 
tiprotons in an electron gas will have completely differ- 
ent convergence behaviors. If this classical example is 
to be representative, protons would then exhibit uniform 
convergence, while antiprotons should exhibit alternating 
convergence. 

Note, finally, that in the problem at hand, response 
theory converges even for arbitrarily large Vq and D + 
S. Since for repulsive spheres the functions f((3V) ap- 
proaches a constant value exponentially, it is even pos- 
sible, in this case, to obtain a quantitative value for the 
limit Vq — > oo by keeping only a finite number of pertur- 
bation terms. 



a uniform, neutral, and interacting jellium 30-32 at a tem- 
perature much lower than the Fermi temperature. This 
hydrogen-in-jcllium problem can be linked to real sys- 
tems in two ways. First, it can be seen as a first step 
in a formal expansion including three- and many-center 
terms which, if carried to all orders, should yield the 
exact total energy of the system, within the adiabatic 
approximation. In can also be used within an effective- 
medium approach 33 in an attempt to take into account 
the many-ion effects in an approximate way. In both 
cases the pair potentials can be used to derive phonon 
spectra. In particular the pair potentials obtained from 
quadratic response were used to predict the infrared and 
Raman vibron frequencies 17 . 

To establish a basis of comparison, we first obtain an 
estimate of the importance of nonlinear corrections to 
the one-atom density, 5p a , using the ab initio program 
VASP 34,35 . This allows us to estimate the effect on the 
pair potentials of the higher-order terms in the determi- 
nation of 5p a . 

To obtain this, we used a cubic cell of side 13.5ao 
containing 74 electrons (for r s = 2), together with a 
6x6x6 fc-space grid and a standard projector augmented 
wave (PAW) pseudopotential for hydrogen 36,37 with cut- 
offs up to 450 eV. The generalized gradient approxima- 
tion (GGA) to the exchange-correlation potential, as pa- 
rameterized by Perdew and Wang 38 , was used, together 
with the Methfessel-Paxton smearing 39 , with a smear- 
ing temperature of a = 0.2 eV. Since it was shown 18 
that the proton-proton pair potential does not depend 
strongly on the choice of a pseudopotential, cutoff en- 
ergy, and exchange-correlation functional, our choice of 
parameters should yield sufficient precision for our pur- 
poses, even considering possible short-range distortions 
due to core overlap. 

The SSPM pair potentials thus obtained are shown in 
Figure 3, and are compared with the VASP results for the 
pair potentials from Bonev and Ashcroft 18 , the quadratic 
response results, and the SSPM pair potentials corrected 
by the inclusion of equation (13) to lowest order in the 
perturbing potential. 

To obtain the latter, we start with (13) and use 

6p^ L (k) ~ I £ X (2) (-k,k',k")^ t (k')K, t (k"). 

k',k" 

This corresponds to the energy diagram shown as Fig- 
ure 5(a). 

We use the second-order noninteracting response func- 
tion of the homogeneous electron gas, Xo i (k,k',k'), in 
the form given by Milchev and Pickenhain 10 . The inter- 
acting response function is then written as 



The hydrogen molecule and connections with -,/( 2 Vk 1/ h-'\ 

n.. u. :i i„ r ,1 u ,.(2)A, l,/ 1,^ _ AO l K ; K > K J 



the Heitler-London approach x s , k',k') = 



e(k)e(k')e(k") : 



We move on to the more realistic system composed of which corresponds to the random phase approximation, 
two initially bound proton-electron systems immersed in This approximation is expected to improve as the density 



10 



SSPM-with d> J 
Ted 

. SSPM-withOLJt 4> red 

1 Ab-initio pair potentia 
Quadratic response 




.* # 



Distance(a Q ) 



FIG. 3: (color online) Comparison of the SSPM results with 
and without contributions from equation (13) to quadratic 
response and ab initio results (from Bonev and Ashcroft 18 ) 
for proton-proton pair potentials at r s = 2. 




Distance(a ) 

FIG. 4: (color online) Contribution of the diagram from Fig- 
ure 5(a) to the proton-proton pair potential for various values 
of r s . 



is increased (see, e.g., Pines and Nozieres 22 ), and should 
be sufficient to give us information on the general be- 
havior of the potential. Finally we choose the Vashishta- 
Singwi form 40 for the local field correction G(k). 

Under these assumptions we obtain the contribution 
to the pair potentials displayed in Figures 3 and 4. 

We see that the inclusion of the potential 4> re d from 
equation (13), and therefore of the density-density in- 
teraction v(q) = e(q) [v c (q) + ^i(q)] , leads to a contri- 
bution to the pair potential that is repulsive at typical 
proton-proton separation. This behavior largely arises 
from the Coulomb repulsion term v c . This contribu- 
tion therefore explains part, but not all, of the discrep- 
ancy observed between the pair potentials obtained from 
quadratic response 17 and ab initio methods 18 . In partic- 
ular, all perturbation-based methods yield a local min- 




(b) 



FIG. 5: (a) Lowest-order diagram contributing to the interac- 
tion (13). (b) Nonreducible diagram equal to diagram (a) in 
the limit R a t — > 0, for identical ions, here protons (diagram- 
matic conventions are explained in Appendix I and Figure 
7 ). 
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FIG. 6: (color online) Contributions of the Coulombic and 
exchange parts to the effective electron-electron potential v, 
for r s = 2. 



imum around R — 1.4ao, whereas no such minimum is 
observed in the VASP pair potential at the considered 
density; such a minimum appears for VASP pair poten- 
tials only at lower densities, corresponding to r s > 3 (see 
Bonev and Ashcroft 18 ). The consequences of the dis- 
appearance of this local minimum at higher density for 
the stability of the hydrogen molecule (or crystal) has 
also been discussed in Nagao et al 17 and Di'ez Muiho and 
Salin 32 . 

Note that for each diagram of the form shown in Figure 
5(a) there are two similar diagrams of the form shown 
in Figure 5(b) which are not reducible and hence not 
contained in (13). In the limiting case where the inter- 
particle distance tends to zero, though, diagrams 5(a) 
and 5(b) should provide the same contribution. There- 
fore diagrams of the form 5(b) should also have an over- 
all repulsive behavior, contributing to further reduce the 
discrepancy between perturbative and VASP results. 

Note also that if we consider separately the contri- 
butions of v c (q) — e(q)v c (q) and /2i = t{q)^i{q), we 
find that the contribution of Jl\ is mostly attractive, as 
can be seen in Figure 6 for r s — 2. We also observe 
that terms arising from equation (13) exhibit very weak 
Friedel oscillations. This distinguishes these terms from 
the other contributions calculated here or ab initio pair 
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potentials 18 (compare Figures 4 and 3), and we conjec- 
ture that this arises from the reducibility of the energy 
diagram. 

We draw the reader's attention to the close similarity 
of the pair potential of equation (15) with corresponding 
terms in the Heitler-London (HL) picture of the isolated 
hydrogen molecule. We can identify in both methods (i) 
the ion-ion repulsions, (ii) the attractive interactions be- 
tween ion b (a, respectively) with the density induced by 
a (6, respectively), (iii) the Coulombic repulsion between 
the one-atom electronic densities, and (iv) an attractive 
exchange contribution from the electrons. The fourth 
term in (15), which has no equivalent in the HL picture, 
goes to zero in the low density, free molecule case. It is 
interesting that two such different approaches, one being 
variational in essence and the other one perturbational, 
yield such similar results. 

Some features of the HL pictures cannot be observed 
in the SSPM, though, because of the nonzero average 
electronic density considered in the SSPM. For example, 
while the pair potential between protons immersed in 
spin-polarized or spin-unpolarized electrons differed only 
by the exchange term in HL, this is no longer the case in 
the SSPM approach. The densities Spi induced by single 
protons are indeed quite different for spin-polarized and 
spin-unpolarized electrons. 



V. CONCLUSION 

To obtain effective interactions between particles im- 
mersed in a well understood system, such as ions im- 
mersed in a uniform jellium, the traditional perturbation 
approach treats all immersed particles as a single per- 
turbation. One finds the free energy of the perturbed 
system as a function of the total perturbation potential. 
The total potential is then separated in the sum of its 
constituents, which allows a natural separation of the 
effective potential in pair-, triplet-, and many-body po- 
tentials. 

Here we have discussed two alternative approaches. In 
the first approach, we start by immersing a single par- 
ticle in the well understood system and then calculate 
the response functions of the new, but perturbed sys- 
tem. Since these response functions are directly related 
to the electronic density or density-density correlation 
functions, they can be obtained by a variety of tech- 
niques, including perturbation theory, simulations, and 
density functional theory and even deduced from exper- 
iment. The remaining particles are then treated as a 
further perturbation of this already perturbed system. 

Such a two-step process is ideal when only one per- 
turbation is strong 24 , but it is intrinsically asymmetrical 
and not ideal when both perturbations are large enough 
to induce nonlinear effects. We have suggested a way 
to improve and symmetrize this procedure which allows 
one to treat an increased number of perturbation terms 
with little additional effort. More specifically, by using 



results of the asymmetrical approach that are exact to 
linear order in the perturbation potential, it is possible 
to construct a symmetrized result exact to quadratic or- 
der and including additional higher-order terms as well. 
We applied this method to two simple noninteracting test 
systems and found that it improves upon standard linear 
and quadratic response, as expected. 

More importantly, it was shown that this simple 
method could be naturally refined by the inclusion of 
higher-order terms describing, in particular, electron- 
electron interactions. These higher-order terms also have 
intuitive physical meaning, and we used this to argue that 
the standard termination of the perturbation series at a 
given order is not the best strategy. The inclusion of the 
higher-order terms was indeed shown to improve con- 
siderably the agreement between the perturbation and 
density functional theory approaches in the problem of 
proton-proton pair potentials. 

Even though we considered here effective interactions 
between identical particles only, the SSP method is par- 
ticularly well suited to the description of systems with 
multiple (say N) types of particles. It indeed reduces the 
computational difficulty from the determination of N 2 /2 
pair potentials to that of finding TV (typically symmet- 
ric) induced densities, from which the pair potentials can 
be obtained in a straightforward manner. It can also be 
generalized to many-center interactions and to magnetic 
perturbations (see Appendix II). 

Finally, a striking similarity is found between terms in 
the pair potentials arising from this approach and from 
the Heitler-London variational approach for diatomic 
molecules, leading to a possible natural generalization of 
the Heitler-London approach to metallic systems. This 
similarity could be explored further by comparing the 
pair potential for a pair of atoms in a jellium, as de- 
rived here, to a pair of atoms in a Wigner-Seitz spherical 
cell ensuring the same average density, which could be 
treated within a Heitler-London-like approach. 

A more detailed comparison of the SSPM and ab ini- 
tio pair potentials for hydrogen and other materials, es- 
pecially at densities higher than those considered here, 
would be a logical next step to this work, as would be a 
more detailed treatment of the many-body interactions, 
both from the ab initio and SSPM perspectives. 
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Appendix I: Derivation of equation (13) 

We first use the Hohcnbcrg-Kohn-Sham approach 41 ' 42 
(and the finite-temperature extension due to Mermin 43 ) 
to write the induced density in terms of the response 
functions of a noninteracting system in a modified ex- 
ternal potential (see, e.g., Lundqvist and March 44 ). We 
write the density as 

«5p(k) = ^x^(-k,k')r(k') 

i k _ (29) 

+ ^E^ 2) (-k,k',k")r ( k')r(k") + ---, 

k',k" 

where Xo™"* are the response functions of the noninteract- 
ing system and T is the effective perturbation potential. 
We write T as 

T(k, A) - XV(k) + c(k, [/>]) - c(k, [po]), (30) 

where /3c(k, [p]) is the first order direct correlation func- 
tion. 

The parameter A is introduced to keep track of the 
order of the expansion. For example, we can write 



where AFi does not depend on the relative positions of 
the perturbation sources and will therefore not contribute 
to the pair potentials. 

Using equations (29), (30), and (32), we can represent 
each term in (34) by a diagram using three types of build- 
ing blocks: namely Xo™^ an d Cq"^ for n > 1, and V ext - A 
few examples are shown in Figure 7. 

Each V ext is connected to a sing le x^.Eachc^ is con- 
nected to n + 1 different x^ ■ Each Xo™^ 1S connected to 
n+ 1 blocks that can be cither V ext , or a Conversely, 
every treelike diagram obeying these rules corresponds to 
a term in equation (34). Note that since Sp(q) obeys the 
same rules, it can be represented by the same diagrams, 
with only one V ex t removed. 

We call reducible those diagrams in (34) that can be 
separated, by the cutting of a single line, into two 
diagrams that depend exclusively on V a or on Vt, respec- 
tively, and that are linear in neither V a nor In Figure 
7, diagram (c) is reducible, while the others are not. 

We want to show that the set of such diagrams is equiv- 
alent to those described by equation (13). First we will 
show that there is a one-to-one correspondence between 
reducible diagrams in (34) and diagrams in (13). Then 
we will show that the prefactors also agree. 



r(k) = Ar«(k) + A 2 r( 2 '(k) + --- 



(31) 



The variation 6c(k, [p]) = c(k, [p]) — c(k, [p ]) can in 
turn be written as a functional expansion: 



5c(k)=^4 1) (-k,k')Mk') 

k' 

+ 4 2) (-k,k',k")^(k')<5p(k") + --- 



k'.k" 



(32) 



Note that here fic^ is the (i + l) th direct correlation 
function of the unperturbed system. This slightly un- 
usual notation is chosen to emphasize the similarity be- 
tween equations (29) and (32). 

Finally the induced density can also be expressed in 
powers of the external potential: 



Sp(k) = XSp^{k) + X 2 5p^(k) + 



(33) 



The variation in the free energy is obtained as before 
through equation (8), which we now write as 



l V P W (k)^t(-k) 
n f-^ n + l 



k,n 



AF ^ + n^ WTi ' 



(34) 



A. Diagrammatic equivalence 

Consider a reducible diagram D contributing to the 
pair potential. Then consider the set S of all diagrams 
that are different from D only by the number of separat- 
ing interaction lines. The separating interaction lines can 
only be connected by Xo^ loops. By summing over the 
different numbers of Xo^ loops and the momenta associ- 
ated with these loops, we can therefore obtain a dressed 
propagator in terms of a dielectric function, 



c o 



(k,k') 



(k,k') 



e(k,k') 



(35) 



k,n 



which takes a simple form in the diagrammatic language 
(see Figure 8). 

All diagrams in S arc therefore included in 
J2{k} ^M{k})A(ki,kj) where Di({k}) is the value, be- 
fore the summation over the momenta {k}, of the dia- 
gram in S with a single separating propagator, k^ 
and kj are the momenta of the separating propagator. 

Since we can construct diagrams contributing to Sp(k) 
with the same building blocks and the same construction 
rules as for diagrams contributing to the pair potential, 
it is easy to find an expression involving only dpi and Cq 1 ^ 
that includes all reducible diagrams contributing to the 
pair potential. An example of such an expression is 
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FIG. 7: Some typical diagrams from the free energy expansion contributing to the pair potential. Loops with n+ 1 legs represent 
response function of order n. The correlation function Cg™' is represented by a solid circle surrounded by n + 1 wavy lines and 
the external potential by solid lines.. Diagram (a) is a quadratic response contribution. Diagrams (b) and (c) represent some 
third-order contributions taken into account by equations (12) and (13), respectively. Diagrams (d), (e), and (f) are examples 
of cubic-response contributions which are not taken into account in the approach presented here. 

+ ••• 



FIG. 8: Diagrammatic expansion for the dressed Sq 1 ' in terms of the bare c|, ij and linear response functions. 



(i) 



k,k' 



Since both (5/9 a (k) and 5pb(k) have an arbitrary num- 



ber of loops on the leg with momenta k, though, 
this expression amounts to screening the Cq 1 ^ interaction 
twice. To take care of this, we can write instead 



(R Q , Rh) 



NLU 



(k,R a )- 



(k,k') 



k,k' 



e(k,k') 



<5^(k',R h ), 



(36) 



where Sp^ LU (k, R a ) is defined as the sum of all diagrams 
in 5p^ L (k,R a ) without \^ loops on the leg with the k 
momentum (U stands for unscreened). To each reducible 
diagram contributing to the pair potential corresponds a 
diagram in (36), and vice- versa. 
For homogeneous systems, we have 



X W(k,k') = x (1) (fc')<S k ,-k', 

4 1) (k,k / ) = 4 1) (fc04,-k', 

i 

1 



and 



e(k,k') e(k) 



(l-c^fe)*^)) 



6p a M(k,R a , b ) =p a , b (k) e 4k ' R -S 
so that (36) simplifies to 



V«*(Ra6) - ^^ i (k) C W(fc)e(fc)^ i (-k) e * R -, 



which is precisely equation (13); we have shown that 
equation (13) contains exactly the reducible diagrams 
from equation (34). Now we need to show that the pref- 
actors of these diagrams also agree. 



B. Prefactors and symmetries 

A diagram D is said to possess a symmetry of order 
m if it is possible to cut m legs from a given xo or c o 
and obtain m identical cut-down parts. Each such sym- 
metry contributes a factor 1/m! to the total prefactor of 
the diagram, as compared to an equivalent asymmetrical 
diagram. 

In order to obtain this "equivalent asymmetrical dia- 
gram" , we replace V ext by V ext — Y^j=i v j i n t ne original 
problem, where n+ 1 is the number of external legs of dia- 
gram D. We then consider the diagram D that is identical 
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to D, apart from its external legs which are all connected 
to different Vj. By construction, diagram D can have no 
symmetry. 

Contributions to diagram D can be obtained, in equa- 
tion (34), by replacing V ex t by any of the Vj and by cal- 
culating the prefactor Pj of the density diagram obtained 
by removing Vj from D. The prefactor Pg of D is there- 
fore 



tron gas can be modeled by the Hamiltonian 45 

where and ~o % are the position and spin operators for 
electron i, while Rj and S 3 are the position and spin of 
the j th magnetic perturbation source. 
The Hellman-Fcynman theorem reads 



n+l 



a external 
■'■ leg 



(37) 



F\ - Fq — j d\{H inti x)\ 

pi r ( 38 ) 

= JJ d\J2j dr/(r - Rj)S l ■ m(r) A , 



Since there is no symmetry in the diagram, it is rela- 
tively straightforward to find its prefactor: each Xo } or whcrc thc local magnetization m is given by 

Cq contributes i\ to the global prefactor. Therefore all 
the Pj are equal, the (n + 1) factors cancel out, and we 
simply get P^ = Pj. By a similar argument, we obtain 



P P ■ — P- P- 



where D a and D b are the separated diagrams. Asymmet- 
ric reducible diagrams therefore have the same prefactor 
in (34) and (13). 

This result can be extended to symmetrical diagrams 
in a straightforward manner: if the two identical dia- 
grams obtained by cutting two legs from a given \ or M 
in a symmetrical reducible diagram, they cannot contain 
the separating interaction line. Therefore all symmetries 
are contained within the separated diagrams. Since each 
symmetry contributes a factor 1/m! to the global prefac- 
tor, to being the number of identical branches involved in 
this symmetry, the symmetry contributions to (13) and 
(34) are the same. 

We have therefore derived the result that there is a 
one-to-one correspondence between reducible diagrams 
contributing to the pair potential [in equation (34)], and 
diagrams contained in equation (13). Since the prefactors 
of these diagrams in each expression also agree, equation 
(13) allows the calculation of the pair potential associated 
with all reducible diagrams. 



Appendix II: Generalizations 

In this appendix we discuss two possible generaliza- 
tions of the methods introduced above: namely, the case 
of magnetic perturbations and that of many-center inter- 
actions. 



A. Magnetic perturbations and the RKKY 
interaction 

The interaction between magnetic perturbations (e.g., 
nuclear spins or magnetic impurity atoms) and an elec- 



m(r) A = ^«y(r,-r)-?^ . 

By analogy with the induced density case, we can de- 
scribe the perturbation using 

V j (r) = jJ2f(r-Ri)S i j . 

i 

Assuming a nonmagnetic unperturbed state, 

m l (r) x = J drxij(r,r')\Vj(r') 

+ J dr'dr" X g(r,r',r")A 2 ^-(r')V fe (r") + --- 

If we have two magnetic perturbing sources (a and b) , 
we find an effective interaction of the form 

0(R a , R 6 , S a , S fc ) = J drV a (r) ■ m 6 (r, R fc , S b ) 

+ V\v) -m a (r,R a ,S a ) -J drdr' Xij (r, r')^ a (r)lf(r'). 

(39) 

Here, as before, V 1 describes the perturbation associ- 
ated with source i and m 4 (r, R i7 S a ) is the magnetization 
induced at r by the presence of a single perturbation of 
spin S z at Ri. 

In the particular case of pointlike magnetic perturba- 
tions [/(r) = <5(r)], the effective interaction between the 
magnetic perturbations takes the simple form 



0(R Q , R b , S a , S 6 ) - JS a ■ m h (R a , R h , S b ) 

+ JS b -m a (R b ,R a ,S a )-J 2 Xt3 (R a ,R b )S?S<>. 

This expression, which goes beyond quadratic response 
[but does not include corrections of the form (13)], could 
be used to improve upon the standard RKKY poten- 
tial, which takes into account only the linear order in 
response 3 . 
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B. Many-center potentials 

Since the sources of the perturbations in equations (12) 
and (13) have not been specified, a straightforward way 
of treating many-particle interactions in this formalism 
is to take either or both sources to be an ensemble of par- 
ticles. This might be especially appropriate for problems 
involving the diffusion of well formed molecules. In the 
hydrogen problem, this could also be used to study the 
molecule-molecule interactions near or beyond the onset 
of metallization. 

Many-center interactions can be treated in a more 



symmetric way through conventional response theory 
(see, e.g., Porter 46 ), but also through a SSPM approach. 
Many-center interactions can also be treated in a more 
symmetric way. The lowest-order TV-body term arises at 
the level of (N — l) th order response. This term is lin- 
ear in the potential from all TV particles. The next most 
dominant contributions will be from terms that are non- 
linear in the potential arising from one particle and linear 
in that arising from the remaining N — 1 . 

We can directly generalize equation (12) to an N-body 
potential treating all terms linear in N — 1 particles ex- 
actly by writing 



N N 

^ )({R}) " (N^m £ £ ^(k)%](-k)e^ - ^ £ X (JV - 1} (k l5 ,k w ) n Vji^e^. (41) 

1 ' 3 = 1 k 1 ' {k,} j=l 



Here Spy] (k) is the density induced by all involved par- 
ticles except j and Spy-i (k) is the component of this den- 
sity that depends on the position of all N — 1 particles 
involved. The latter definition is required simply to avoid 
double counting the potentials involving less than N par- 
ticles. 



Equation (41) provides in principle a way to obtain an 
expression that is exact up to N th order, requiring only 
the explicit knowledge of the (N — \) th order response 
function. On the other hand, it requires knowledge of the 
density induced by N groups of iV — 1 particles, which 
rapidly becomes more difficult when N > 2. 
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